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摘 ”要 : 基 十 完全 重合 型 区 域 分 解 技 巧 ， 本 文 提出 了 一 种 求解 非 定常 Stokes 方 程 的 有 限 元 并 行 算 法 。 该 
算法 的 基本 思想 是 首先 对 空间 施行 完全 重 玖 型 区 域 分 解 ， 然 后 各 个 处 理 器 使 用 向 后 Euler 格式 独 
立 并 行 求解 关 填 时 间 t 的 常 微分 方程 ， 在 整个 关于 时 间 的 迭代 过 程 中 ， 无 需 处 理 器 间 的 通信 ， 具 
有 良好 的 并 行 性 能 。 该 算法 中 每 个 处 理 器 所 负责 的 子 问题 是 一 个 全 局 问题 ， 它 定义 在 整个 求解 区 
域 上 ， 但 绝 大 部 分 自由 度 来 自 其 所 负责 的 了 区域 ， 从 而 使 得 该 算法 稍 加 修改 现 有 的 串 行程 序 即 可 
实现 相应 的 并 行 计算 ， 实 现 简单 ， 具 有 重要 的 使 用 价值 。 同 时 通过 数值 算 例 ， 在 曙光 集群 并 行 机 
上 编程 实现 了 上述 算法 ， 验 证 了 其 有 效 性 。 
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RN ÆR? (d = 2,3) 中 具有 Lipschitz 连 续 边 界 的 有 界 开 集 ， 我 们 考虑 非 定常 Stokes 方程 


ue — vAu +Vp= f, €Q x (0,T] 内 ， 


div u = 0, t Q x (0,T] A x 
u=0, 在 aN x (0,T] E, 
u = Uo; EQ x {0} 上 ， 


Hp u = (u(x,t), uala, t) 为 速度 向 量 ，P = pla, t) AEH, f = (filz, t) , fala,t)) X 
ARI, v 是 粘性 系数 ，wo 是 初始 速度 向 量 ， 开 是 一 给 定 的 有 限时 间 ，w = St. J0 是 流 
体力 学 的 基本 方程 一 Navier-Stokes 方程 的 线性 化 形式 ， 描 述 了 蠕动 流 的 流动 情况 。 因 此 ， 其 并 
行 数值 求解 方法 的 研究 具有 重要 的 现实 意义 。 

区 域 分 解 算法 是 一 类 数值 求解 偏 微 分 方程 的 有 效 方法 ， 它 基于 “分 而 治之 ”的 思想 ， 将 求 
解 区 域 分 解 成 若干 个 形状 相对 较 规则 的 子 区 域 ， 使 得 原 方程 转化 为 在 子 区 域 上 求解 ， 从 而 将 大 
型 问题 分 解 为 小 型 问题 、 复 杂 边 值 问 题 分 解 成 简单 边 值 问 题 、 串 行 问题 分 解 成 并 行 问题 ， 以 达 
到 扩大 求解 问题 的 规模 、 提 高 计算 精度 、 加 快 计算 速度 的 目的 四。 对 于 非 定 常 问题 ， 传 统 的 区 
域 分 解 算法 是 首先 对 时 间 进 行 离 散 ， 然 后 在 每 一 时 间 步 对 空间 施行 区 域 分 解 实 现 并 行 计算 。 这 
”收藏 日 期 : 2009-12-01。 作者 简介 ， 尚 月 强 (1976 年 12 月 后), 男 , 博士， 副教授 . 研究 方向 ， 偏 微分 方程 数值 解 与 

并 行 计算 . 
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类 方法 的 一 个 最 大 缺点 是 处 理 器 之 间 的 通信 较 频繁 ， 因 而 其 并 行 性 能 往往 不 够 理想 。 本 文 提出 
求解 非 定常 Stokes 方 程 (1) 的 一 种 新 型 区 域 分 解 有 限 元 并 行 算法 ， 该 算法 首先 对 空间 施行 完 
重 玲 型 区 域 分 解 ， 然 后 在 每 个 处 理 器 使 用 向 后 Euler 格式 独立 并 行 求解 关于 时 间 t 的 常 微分 方 
程 ， 在 整个 关于 时 间 的 迭代 过 程 中 ， 无 需 处 理 器 间 的 通信 ， 具 有 良好 的 并 行 性 能 。 同 时 我 们 在 
曙光 集群 并 行 机 上 进行 了 数值 试验 ， 验 证 了 其 有 效 性 。 


2 ”预备 知识 


2.1 非 定常 Stokes 方 程 解 的 存在 了 崔 一 性 

对 于 非 负 整数 上 ， 我 们 按 通常 意义 定义 Sobolev 空间 五 R(Q)Q， 并 令 丽 (Q) = {u € HHN): 
ulaz = 0}， 这 里 wleo = 0 是 在 迹 的 意义 下 。 本 文中 ， 我 们 用 (-,-) KL) (d = 1,2,3) 上 的 
标准 内 积 ， 用 jp 表 五 &(O)( > m) 空间 中 的 m- 全 范 。 令 


X = HQ), Y =Q), M = LQ) = fa € LQ): J qdz = 0 \, 
Q 
并 定义 双 线 性 形式 a(., .), ad(., :) 
a(u,v) =v(Vu, Vv), d(v,q) = (divv,g), Vuve x, q€M, (2) 
我 们 可 得 方程 (1) 的 变 分 形式 为 : 求 (w,p) EX x M， 使 得 对 任意 给 定 的 te (0,7), A 
(uz, v) + a(u,v) — d(v,p) + d(u,q) =(f,v), VvV(v,9)EXxM. 
u(0) = uo. 


假设 初 值 uo 及 右 端 项 了 满足 以 下 条 件 : 
HO up € (H?(Q) HG (M))*, divuo =0, f, fe € L?(0, T; L7(Q)*) 满足 


T 1/2 
luole+( f (A+B) <e (4) 


本 文中 c 是 一 大 于 零 的 常数 ， 它 在 不 同 的 地 方 可 能 代表 不 同 的 值 。 对 于 方程 (3)， 我 们 有 以 下 结 
果 。 

E2104 ”假设 HO 成 立 ， OH RIPRA CH (k > 1) 型 边界 的 有 界 区 域 或 Ri 中 的 凸 
多 边 形 或 凸 多 面体 区 域 (k = 1)， 则 对 于 任意 的 0 <T < co， 方程 (3) 存在 唯一 解 (w,p) 满足 


Usup (eO + PON + fet) < 7 
T 
sup _o(t)||uell} +f a(t) ( luill? + llpell? + lunlia )dt < c, (6) 
O<t<T i 


Ht o(t) = min{1,t}. 

2.2 ETE Stokes FAVA PR ITIBIL 

ETN) = {KENKA — ARIE IEA EET) 的 尺度 为 产 的 网 格 前 分 ，X8(O) C 
X, MRQ) c MAAN FART (Q) 的 两 个 有 限 元 空间 ， 满 足 蜗 近 性 质 : 任 给 (w,p) € 
H+A(Q)4 x HE(Q)(k > 1)， 存 在 (rnuw pnp) € XP(Q) x M2(Q)， 使 得 


Ju — mhuli < ch*}julligs.o, lp—papllon < chsllplls,no, OSs <k, (7) 
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和 inf-sup (LBB) 条 件 ， 存在 8 > 0， 使 得 


(div v, q) 
sup ; 
vex? (Q), w0 ||Yvllo,n 


Bllgallo,n < vg € M(N). (8) 


使 用 向 后 Euler 格式 对 时 间 进 行 离散 ， 我 们 可 得 方程 (3) 的 全 离散 有 限 元 格式 为 ， 求 (w?,p?) 
E XR(Q) x MR(Q9)， 使 得 


Ay (uk — up v) + a v) ~ d(v, pR) + dlug, a) = (f",v), Y (v,a) € XRO) x MP(Q), 
u? = Uho; 


其 中 At 为 时 间 步 长 满足 0 < At <1, P= hfe, Jdt, tp =nAt, n=1,2 N, N= 
[T/At], uno 是 wo WAA RER. 

对 于 方程 (9)， 我 们 有 以 下 结果 。 

定理 2.284 在 定理 2.1 的 条 件 下 ， 对 所 有 的 m =1,2,---,N, 方程 (9) 存在 唯一 解 (wu?,p%) 
€ XP(Q) x MP(Q) 满足 


IV (ulta) — up )ilo,0 + lp(tn) 一 焉 oo < e(At +h), 1<s<h, (10) 


其 中 一 c(a(t), Q, uo, f, T)» 


3 ARH Stokes 方程 基于 完全 重 到 型 区 域 分 解 的 有 限 元 并 行 算法 


3.1 ”完全 重合 型 区 域 分 解 技 巧 

首先 将 求解 区 域 人 分解 成 互 不 重合 的 子 区域 D1, D2,… , DJ， 然后 在 人 9 内 将 每 一 子 区 域 向 
外 扩展 一 定 尺 度 获 得 mi,Q2,…. ,QJ (Q cR, J= 1,2,--- JJ), 这 些 子 区 域 Q ,92 …， 7 构 
成 Q 的 一 个 重 基 型 区 域 分 解 。 每 台 处 理 器 负责 一 个 子 区 域 ， 各 自 对 所 负责 的 子 区 域 生成 
尺度 为 h 的 细 网 格 Th(Q;)， 并 对 其 余 区 域 生成 尺度 为 万 > h 的 粗 网 格 TH(Q\Q;)， 然 后 
采用 自 适应 处 理 等 技巧 使 得 粗细 网 格 在 交界 处 相 容 ( 即 任 两 个 单元 的 交集 或 者 为 空 集 ， 
或 者 为 一 公共 顶点 ， 或 者 为 一 公共 边 或 为 一 公共 面 )， 记 这 样 的 整个 区 域 的 多 人 尺度 网 格 
Hy TZ" (Q) (7 = 1,2,… ,J)。 这 样 的 网 格 也 可 通过 对 整个 区 域 的 初始 粗 网 格 TH(Q) 施行 局 部 
加 密 并 通过 自 适应 过 程 使 之 相 容 得 到 。 图 1 描绘 了 四 个 子 区 域 情形 的 完全 重 敬 区 域 分 解 和 相应 
的 网 格 剖 分 情况 。 在 该 技巧 中 ， 相 应 于 每 一 子 区 域 07 的 子 问题 定义 在 整个 求解 区 域 人 上， 但 
HAMS ARK EAP RIM; 每 一 个 子 问题 事实 上 是 一 个 全 局 问题 ， 因 而 稍 加 修改 现 
有 的 串 行程 序 即 可 实现 并 行 计 算 ， 实 现 简单 ， 通 信和 需求 少 。 

3.2 ”基于 完全 重 佛 型 区 域 分 解 的 有 限 元 并 行 算法 

设 Di, D2,… , DJ 是 求解 区 域 Q 的 一 个 互 不 重 登 的 区 域 分 解 ，0 CQ = 1,2,… ,J) 通 
过 扩展 Dj (7 = 1,2,… ,J) A. TOG = 1,2,… ,J) 是 的 多 尺度 网 格 剖 分 ， 它 在 Qj 内 
的 尺度 为 h， 而 在 远离 Qj; 的 地 方 尺度 为 五 (参见 图 1)， 这 些 T8"(Q) (7 = 1,2,… ,四 构成 
了 的 一 个 完全 重 辣 型 区 域 分 解 。 设 XX; (O) 和 M9; jy(Q) EAT TOO (7 = 1,2,… J) 
的 有 限 元 空间 ， 则 求解 方程 (1) 基于 上 述 完 全 重生 型 区 域 分 解 的 有 限 元 并 行 算法 如 下 : 

算法 1 非 定 常 Stokes 方 程 基于 完全 重 秋 型 区 域 分 解 的 全 离散 有 限 元 并 行 算法 

对 n= 1,2,.……,N;: 
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1) HAE (wh, PH) < Xal) x ME ng (O) (i = 1,2,… ,J)， 使 得 
ub” eget , K n 
(E v) + alui v) — dopin) + dukia a) = (F° 0), 


At 
v (v,q) € X Pr ng (Q) x Mfr n (9), 
AA = Pyh juo. 


2) E Dj 内 取 (uh, ph) = (ui ys Pin) (j = 1,2,.… ,J)。 其 中 Pung :Yo X Fr n,j(Q) # JE 
交 投 影 ， 定 义 为 


(PH hju v) = (u,v), VuEY, ve Xha), j=1,2, ,J. 


A 


X 


VA/ 


Š 


RERA N 
ANNS 
wa 


LVN 


oO 
N 


X 
A 


Bl: ”四 个 子 区 域 情形 完全 重 登 型 区 域 分 解 


3.3 ”算法 分 析 

由 算法 1 可 知 ， 新 算法 中 每 一 子 问题 实际 上 是 一 个 全 局 问题 ， 从 而 稍 加 修改 现 有 的 串 行 程 
序 即 可 实现 并 行 计算 ， 实 现 简单 ， 在 关于 时 间 的 和 迭代 过 程 中 ， 无 需 处 理 机 之 问 的 通信 ， 有 具有 较 
好 的 并 行 性 能 。 定 义 分 片 范 数 


J 1/2 J 1/2 
Nv = (Dv,) > Ie- loe= (X le-o) - 
j=1 j=1 


针对 定常 Stokes 方 程 ， 基 于 两 重 网 格 离散 技巧 ， 何 银 年 教授 等 人 在 文献 [8] 中 提出 了 类 似 于 完 
全 重 秋 型 区 域 分 解 的 有 限 元 并 行 算法 ， 并 在 网 格 和 有 限 元 空间 的 一 些 假 设 条 件 ( 即 ; 逼近 性 、 
反 估计 、 超 收敛 性 、LBB 条 件 ) 下 ， 得 到 下 列 误差 估计 


lv -woa + Ilp = p” lllo < c(h? + H°*"), 1<s<k. (11) 


对 于 非 定常 问题 ， 由 于 与 时 间 相关 ， 其 理论 分 析 比 定常 问题 要 复杂 得 多 。 本 文 我 们 就 所 提出 的 
新 算法 1， 主 要 从 数值 模拟 的 角度 考察 其 有 效 性 及 并 行 性 能 ， 其 相关 理论 分 析 我 们 将 在 后 续 工 
作 中 给 出 。 


4 ”数值 试验 
本 节 我 们 给 出 两 个 数值 算 例 以 验证 算法 的 有 效 性 。 数 值 试验 平台 为 西安 交通 大 学 理学 院 智 


能 信息 处 理 与 计算 实验 室 的 曙光 集群 并 行 机 ， 其 每 个 结 点 具有 8 核 2.0GHz CPU, 2GB x 8 内 
存 ， 各 结 点 通过 20Gbps InfiniBand 连 在 一 起 。 消 息 传递 软件 为 MPICH。 
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算 例 1 解析 解 
针对 求解 区 域 2 = [0,1] x [0,1] 和 已 知 解析 解 的 情形 ， 即 准确 解 为 


uy = £?(x — 1)*y(y — 1)(2y — 1) cos(t), 
ug = —y*(y — 1)?2x(x — 1)(2x — 1) cos(t), (12) 
p = (3x? + 3y? — 2) cos(t), 


FONT BCT BARS TF: PT ACG RR A DR 9) MD SF Re, PF Pb — Py 元 和 文献 [9] 中 
所 讨论 的 稳定 化 已 — Py F Py — Po Teil SARA PRCA, RIL KT h TREE; RIE 
察 算法 1 的 并 行 性 能 。 

Pib 一 PP 元 的 有 限 元 空间 定义 为 


Xing (Q) = {v € Hy (Q)? : vle € (Pi @ span{b(K)})?, VK Ee TEQ}, j=1,.…,), 
Ming (Q) = {q € LEN) NCN): qk € Pi, YK Ee TEQ) }, j=1,.…,), 


SEP OK) 为 单元 天 EKAR 3 ARA, WEK) =0, P 为 线性 多 项 式 函 数 空间 。 
对 于 稳定 化 的 已 — P CAP, — Po vt, ， 其 速度 空间 定义 为 


X9 ni (Q) = {v E€ HQP :vk € (P, VK € TQR j=l, od: 
压力 空间 分 别 为 
Miz ng (Q) = {4 € LO) NCN): dle E Pi, VK ETP}, j=1, yd, 


Mf n = {a € L6(0) : alk € Po, VK E TQ) }, j = Tiss J, 


相应 的 稳定 化 格式 为 
(BASHA) + aluin o) — dv, phi) + alee a) + OOM.) = 2.0) 
V (v,a) € XE ns(0) x MEn, j (9), (13) 
utp, = PH, juo. 
其 中 稳定 项 C(.，) 定义 为 
C(p,a) = (T - D)p, (I - 1)4), Yp, q € L?’(9), (14) 


H: L?(Q) > Po (Xt Py — Py 7G) RIE: L) > P (Xt Py — PC) ARARBAT. 
FAC» FRAT PRS SC AF DK So} RD AS OS a a AF SR 


Dı = (0,1/2) x (0,1/2), D2 = (1/2, 1) x (0, 1/2), 
Ds = (0, 1/2) x (1/2,1), Da= (1/2, 1) x (1/2, 1), 


SR a CEO. AKG RE — DK BR Jey Ob HER BE OY DK RRO; (0 = 1,2,3,4)， 并 分 别 就 h = 
n-?, H=Vh(n=4,6,8, 10) 47i, Hhyv=1, At=0.001. WHA RWIS. 
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表 1: 算法 近似 解 的 误差 :At = 0.001， 四 子 区 域 情形 


Er | We rile 
AANA Italo 
1/64 | 1/8 0.99330 1.00678 
1/100 | 1/10 1.02712 1.05341 
1/6 1.07247 | 1.32402 
/10 ; 
/4 


Uni 收敛 率 PL2 收敛 率 


Keng ee 
~~ 
oje] 
QIO 


T =0.01 | P -P 


1 0.94602 1.09591 
1.05807 1.00622 
1.01597 1.00692 
0.92408 1.01468 
0.948467 | 1.09713 
0.949655 | 1.08312 
0.990862 | 1.13267 
1/6 1.07118 1.29968 
1/10 0.90981 | 1.10747 


1 


1/16 
1/36 
1/64 
/100 
1/16 


P, — Po 


1/36 
1/64 
1/100 


1/6 
1/8 
1/10 


Pib- Py 


1/16 | 1/4 
1/36 | 1/6 
1/64 | 1/8 
1/100 | 1/10 


P,- Pı 


1/16 
1/36 
1/64 
1 


| 
ll 
> 


1/10 0.87750 1.01567 


计算 结果 表明 : KT = 0.01 时 ,三 种 有 限 元 方法 所 得 近似 解 关 于 网 格 参数 hh 都 取得 了 接 
近 标 准 有 限 元 方法 的 一 阶 收敛 率 ， 但 随 着 时 间 向 前 推进 ， 其 速度 的 精度 和 收敛 阶 有 所 降低 。 
在 所 得 近似 解 精 度 方面 ， 不 论 是 速度 还 是 压力 ，Pib — Py 元 所 得 近似 解 的 精度 明显 高 于 稳定 
化 的 号 — Pi AP, -西元 ,这 有 别 于 标准 的 全 局 算法 (参见 文献 [9] MAHAR): 而 稳定 化 
的 已 一 请 元 和 户 一 己 元 之 间 ， 记 一 忆 元 所 计算 的 压力 精度 显著 高 于 只 一 局 元 ,但 其 所 得 
速度 的 精度 却 稍 差 于 PP — Po Jie 

其 次 ， 我 们 将 求解 区 域 分 别 划分 成 2-16 个 子 区 域 进行 相应 的 计算 ， 以 考察 算法 的 并 行 性 
能 。 表 2 及 图 2 分 别 报告 了 三 种 不 同 的 有 限 元 方法 在 曙光 集群 并 行 机 上 的 运行 时 间 、 加 速 比 、 
并 行 效率 和 加 速 比 随处 理 器 数 的 变化 情况 ， 其 中 Ves 为 网 格 TP*(Q) (7 = 1,2,… , J) 中 三 角 单 
元 顶点 数 的 最 小 者 ，h = 1/100, H= 1/10, At = 0.001。 数 值 结果 表明 : 算法 取得 了 较 好 的 并 
行 性 能 ， 而 且 关于 时 间 的 迭代 次 数 越 多 ， 并 行 性 能 越 好 ; 当 人 = 0.1 时 ， 记 6b 一 只 元 在 处 理 器 
个 数 为 4、6 时 还 取得 了 超 线性 的 加 速 比 和 并 行 效率 。 
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表 2: 并 行程 序 的 运行 时 间 代 (J) ( 秒 )、 加 速 比 Sp = FQ 和 并 行 效率 Ep = PAB 


子 区 域 数 J | Ves 


7234 | 176.55 138.52 135.56 | 1.0 
se 


2.75 | 0.92 0.92 


fare far af af 
1.02 | 470.58 | 2.88 s 464.63 0.95 
e 1.01 | 346.59 | 3.91 364.01 | 3.65 | 0.91 


2202 | 323.53 252.3 0.88 
1830 | 243.57 ; . 
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(a) T = 0.01 (b) T =0.1 
图 2 h= 1/100, H=1/10, At = 0.001 时 算法 的 加 速 比 FD AGA we Bae ay RACH 


算 例 2 后 台阶 问题 

本 算 例 考虑 后 台阶 间 题 ， 它 是 检验 一 个 算法 有 效 性 的 很 好 的 基准 问题 之 一 ， 其 求解 区 
域 及 边界 条 件 如 图 3 所 示 ， 其 中 体积 力 = 0。 我 们 将 求解 区 域 分 解 成 五 个 互 不 重 亚 的 子 区 
域 ( 见 图 4)， 然 后 将 每 一 子 区 域 向 外 扩展 尺度 为 五 的 区 域 以 获得 相互 重 登 的 区 域 分 解 ， 其 
Ph = 1/8, H = 1/4， 所 使 用 的 有 限 元 为 二 阶 Taylor-Hood 元 。 时 间 步 长 At = 0.01， 初 始 速 
EuL =0, v=1. 

图 5 描绘 了 流体 达到 稳定 状态 ， 即 满足 迭代 终止 条 件 


ubr! — yl 
Nin — Uralon P 210, JA (15) 
lugn log 


时 速度 和 压力 的 等 值 线 。 需 注意 的 是 图 5 中 每 一 子 图 在 z = 0,4,8,12 处 的 垂直 线 是 子 区 
Dj (ij = 1,… ,5) 的 人 工 边 界线 。 该 算 例 进一步 验证 了 本 文 所 提出 的 并 行 算 法 的 有 效 性 。 


(-4. 2) (16, 2) 
ui = O,ug = 0 Ou 
Ag “pty =O) wo 
(-4. 0) uy = 0, u2 = 0 w= 
(0. -1) (16, -1) 
图 3: 后 台阶 问题 
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图 4: 后 台阶 问题 的 区 域 分 解 


(c) 压力 p 


图 5: 台阶 问题 达到 稳定 状态 (了 = 0.92) 时 近似 解 的 等 值 线 


5 结论 


基于 完全 重 伙 型 区 域 分 解 技 巧 ， 本 文 提出 了 一 种 新 的 求解 非 定常 Stokes 方 程 的 有 限 元 并 行 
算法 。 该 算法 实现 简单 ， 可 充分 利用 现 有 的 串 行 软件 进行 并 行 编程 ， 同 时 具有 良好 的 并 行 性 
能 。 数 值 算 例 验 证 了 其 有 效 性 。 
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A Parallel Finite Element Algorithm Based on Fully Overlapping 
Domain Decomposition for the Time-dependent Stokes Equations 


SHANG Yue-qiang’”, HE Yin-nian? 


(1- School of Mathematics and Computer Science, Guizhou Normal University, Guiyang 550001; 
2- Faculty of Science, Xi’an Jiaotong University, Xi’an 710049) 


Abstract: Based on the fully overlapping domain decomposition, a parallel finite element algorithm for 
the time-dependent Stokes equations is proposed. The basic idea of the algorithm is to first discretize 
the spatial space by using the fully overlapping domain decomposition technique, then independently 
solve a system of ordinary differential equations with respect to time by the backward Euler scheme 
in overlapped subdomains. During the time-iterations, there is no communication between processors 
and hence a good parallel performance can be obtained. In this algorithm, each subproblem is a global 
problem with the vast majority of the degrees of freedom associated with the particular subdomain that 
it is responsible for, and hence can be solved with other subproblems in parallel by using an existing 
sequential solver without extensive recoding. Numerical tests on a Dawning parallel cluster illustrated 
the efficiency of the algorithm. 

Keywords: Stokes equations; finite element method; overlapping domain decomposition; parallel 


algorithm 
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